lung ligand-receptor analysis

In this note we exemplify the Ligand-Receptor analysis as was shown in (Cohen et al. Cell 2018).

For further questions please contact Amir Giladi or Dikla Glebard Solodkin

Ligand-Receptor analysis

loading Metacell library

In order to run the following commands you will need mat and mc objects in your database directory and update the following chunk with the correct paths and id’s

library("metacell")
if(!dir.exists("saved_work")) dir.create("saved_work/")
scdb_init("saved_work/", force_reinit=T)

## initializing scdb to saved_work/

tgconfig::override_params("annotations/lung_params.yaml","metacell")
src_functions = "scripts/mc_modeling_functions.R"
if(!dir.exists("results")) dir.create("results/")
scfigs_init("results/")
if(!dir.exists("results/figure3")) dir.create("results/figure3")
if(!dir.exists("results/figure3/pairs")) dir.create("results/figure3/pairs")
if(!dir.exists("results/figure3/pops")) dir.create("results/figure3/pops")
if(!dir.exists("results/figureS3")) dir.create("results/figureS3")
if(!dir.exists("results/figureS3/lr_kinetics")) dir.create("results/figureS3/lr_kinetics")
id = "lung_kinetics_sorted" # here you can update for the name of your MetaCell object
mat_id = id
mc_id = id

Interaction map

Building and projecting an interaction graph of ligand-receptor pairs

#' @param mat_id 
#' @param mc_id
#' @param dir_nm
#' @param fig_fn file name for figure
#' @param frac_threshold choosing only genes which are expressed in more than 15% of the cells of the maximal metacell  
#' @param K number of gene modules (for visualiztion)
#' @param modules_corr_threshold numeric value between 0 to 1, increasing the number affects the density of the graph
#' @param genes_corr_threshold numeric value between 0 to 1, increasing the number affects the density of the graph in each module
#' @param jitter_x random range for x axis
#' @param jitter_y random range for y axis

library(reshape2)
library(Rgraphviz)
library(plyr)
lig_rec_map <- function(mat_id, mc_id,lig_rec_fn ='annotations/ligand_receptor_mouse.csv' ,fig_nm = "results/figure3/interaction_map.png",frac_threshold = 0.1,
                        K=15,modules_corr_threshold = 0.55, genes_corr_threshold = 0.3,jitter_x = 35,jitter_y = 35){
  source(src_functions)
  # load objects
  mat = scdb_mat(mat_id)
  mc = scdb_mc(mc_id)
  color_key = unique(mc@color_key[,c("group","color")])
  name2color = color_key$color
  names(name2color) = color_key$group
  color2name = color_key$group
  names(color2name) = color_key$color
  if(file.exists(paste0(.scdb_base,mat_id,".dus.Rda"))){
    # loading a downsampled matrix (dus) which was already computed
    load(paste0(.scdb_base,mat_id,".dus.Rda"))
  }
  else{
    # creating a gene set from the rownames of mat@mat in order to generate a downsampled matrix (dus)
    all_gene_sets = rep(1,length(rownames(mat@mat)))
    names(all_gene_sets) = rownames(mat@mat)
    scdb_add_gset("all_mat_genes",gset_new_gset(all_gene_sets,"all_mat_genes"))
    set.seed(1436)
    dus = as.matrix(gset_get_feat_mat(gset_id = "all_mat_genes",mat_id = mat_id,downsamp = T))
    # removing rows (genes) with zero umi count
    dus = dus[which(rowSums(dus) !=0),]
    save(file=paste0(.scdb_base,mat_id,".dus.Rda"),dus)
  }
  set.seed(27)
  # removing columns (cells) which are not in the mc object 
  dus = dus[,intersect(names(mc@mc),colnames(dus))]
  # loading ligand receptor table
  proper=function(x) paste0(toupper(substr(x, 1, 1)), tolower(substring(x, 2)))
  rec_lig = read.delim(lig_rec_fn, sep = ",", row.names = 1)
  ligand_receptor = as.matrix(table(rec_lig$ligand,rec_lig$receptor))
  rownames(ligand_receptor) = unlist(lapply(rownames(ligand_receptor),function(x) proper(x)))
  colnames(ligand_receptor) = unlist(lapply(colnames(ligand_receptor),function(x) proper(x)))
  lfp = log2(mc@mc_fp)
  sigs = intersect(intersect(rownames(lfp),rownames(dus)), union(rownames(ligand_receptor), colnames(ligand_receptor))) # list of relevant ligand receptor genes
  # filter out low abundance genes
  mc_mc = mc@mc[intersect(names(mc@mc),colnames(dus))]
  # m is the matrix of the normalized average of UMI counts for each LR gene (rows) in each metacell (columns) - created from the downsampled UMI count matrix, averaged by metaclles size
  m = t(apply(dus[sigs,], 1, tapply, mc_mc, sum)) 
  sizes = table(mc_mc)
  m = sweep(m,2,as.vector(sizes),"/")*min(sizes) 
  # define a "maximal metacell" for each gene and calculating in how many cells of this metacell, the gene is expressed 
  max_m = apply(m,1,max)
  arg_max_m =  apply(m,1,which.max)
  count_m = sapply(sigs, function(x) (sum(dus[x,names(mc_mc[mc_mc == arg_max_m[x]])] >0)/sizes[arg_max_m[x]]))
  png("results/max_avg_m_vs_count_m_ds.png")
  plot(log2(max_m),count_m,col = mc@colors[arg_max_m],xlab="maximal normalized averaged log2 umi count of a gene in a metacell",ylab = "fraction of cells with at least one molecule of the gene in the maximal metacell")
  abline(h = frac_threshold)
  dev.off()
  names(count_m) = sigs
  if(file.exists("saved_work/bad_marks.Rda")){
    load("saved_work/bad_marks.Rda")
  } 
  else{
    bad_marks = c()
  }
  # filter out genes that are expressed in only less than frac_threshold of the cells in their "maximal metacell"
  genes = setdiff(names(which(count_m >  frac_threshold)),bad_marks) # filtered list of relevant ligand receptor genes
  # creating a data frame of ligand-recptor pairs from the filtered gene list 
  ligands = intersect(genes, rownames(ligand_receptor)); receptors = setdiff(genes, ligands)
  rec_lig_2 = melt(ligand_receptor)
  colnames(rec_lig_2) = c("ligand", "receptor", "interaction")
  rec_lig_2[,1] = as.vector(rec_lig_2[,1]); rec_lig_2[,2] = as.vector(rec_lig_2[,2])
  rec_lig_3 = rec_lig_2[ rec_lig_2[,1] %in% ligands & rec_lig_2[,2] %in% receptors & rec_lig_2[,3] == 1,]
  # choosing highly variable genes for each metacell
  nms = choose_genes_from_mc(mc,mat, nms_per_mc=15,bad_genes = c(bad_marks,setdiff(rownames(mc@mc_fp),rownames(dus))),nms_thresh = 2)
  genes = union(ligands, receptors)
  all_genes = union(nms, genes)
  gene_set = rep(1,length(all_genes))
  names(gene_set) = all_genes
  # calculting correlation between genes (list contains ligands, receptors and other significant genes) according to the downsampled umi counts matrix
  # running hierarchical clustring on this correlation matrix 
  feat = dus
  k_nonz_exp = tgconfig::get_param("scm_k_nonz_exp",package = "metacell")
  feat = log2(1+k_nonz_exp*as.matrix(feat))
  feat_cor = tgstat::tgs_cor(t(as.matrix(feat)[all_genes,]),spearman = TRUE)
  hc = hclust(as.dist(1-feat_cor), "ward.D2")
  ct = cutree(hc, K)
  C = feat_cor
  m = apply(as.matrix(feat)[all_genes,], 2, tapply, ct, sum)
  corr_modules = cor(t(log(1 + m)), method = "spearman")
  C3 = corr_modules
  C3[ C3 < modules_corr_threshold] = 0
  N = nrow(C3)
  # drawing the interactions graph, N nodes for each genes module
  rEG <- new("graphNEL", nodes=as.character(1:N), edgemode="undirected")
  e = which(C3 > 0) # edges between similar modules 
  n1 = ceiling((e)/N) # node1 in the edge
  n2 = 1+((e-1) %% N) # node2 in the edge
  # drawing edges between modules
  rEG = addEdge(as.character(n1[n1!=n2]), as.character(n2[n1!=n2]), rEG, C3[e[n1!=n2]])
  g = layoutGraph(rEG, layoutType="neato")
  x_cl = nodeRenderInfo(g)$nodeX # getting x position of each module
  y_cl = nodeRenderInfo(g)$nodeY # getting y position of each module
  names(x_cl) = rownames(C3)
  names(y_cl) = rownames(C3)
  # right now we have positions for the centers of each genes module, now we would like to draw each gene  
  C2 = C; C2[ C2 < genes_corr_threshold] = 0
  C2 = C2 / rowSums(C2)
  # choosing positions for each gene according to similarity 
  x = apply(C2, 1, function(z) sum(z[z > 0 ] * x_cl[ ct[ names(which(z > 0))]])) + runif(nrow(C2),-jitter_x,jitter_x)
  y = apply(C2, 1, function(z) sum(z[z > 0] * y_cl[ ct[ names(which(z > 0))]])) + runif(nrow(C2),-jitter_y,jitter_y)
  gene_coords = cbind(x,y)
  rownames(gene_coords) = all_genes
  IM = lfp[all_genes, ]
  wmax = rep(NA, length(all_genes)); names(wmax) = all_genes
  # for each gene, chossing the cell type with its maxmial expression levels
  mean_lfp = tapply(1:length(mc@colors),color2name[mc@colors],function(x) {
    if(is.vector(IM[,as.integer(x)])){return(list(IM[,as.integer(x)]))}
    else{ return(list(rowMeans(IM[,as.integer(x)])))}})
  mean_lfp = as.data.frame(mean_lfp)
  wmax[rownames(IM)] = colnames(mean_lfp)[apply(mean_lfp, 1, which.max)]
  # display only ligands and receptors genes
  gene_coords = gene_coords[genes,]
  wmax = wmax[genes]
  # plotting the interaction map
  png(fig_nm, height = 2000, width = 2000)
  par(mar=c(5,20,5,0),xpd=TRUE)
  plot(gene_coords[,1], gene_coords[,2], type = "n", axes = F, xlab = "", ylab = "")
  # plot segments (edges) between pairs of ligand-receptor with interaction
  with(rec_lig_3[ rec_lig_3$ligand %in% genes & rec_lig_3$receptor %in% genes,], 
       segments(gene_coords[ligand,1], gene_coords[ligand,2], gene_coords[receptor,1], gene_coords[receptor,2], lwd = 1.5,col = "gray60")) 
  # plot points for each gene, color them by the cell type of the identified metacell
  points(gene_coords[,1], gene_coords[,2], pch = 21, cex = 4, lwd = 4,
         bg  = ifelse(rownames(gene_coords) %in% ligands, name2color[wmax], "white"), 
         col = ifelse(rownames(gene_coords) %in% ligands, "black", name2color[wmax]))
  legend(min(gene_coords[,1])-55,max(gene_coords[,2]) + 15,legend = c("ligand","receptor"), pch=21,pt.lwd = 1,cex = 3,pt.bg=c("gray","white"),box.lwd = 0,bg=NA,pt.cex = 4,title = "LR type",title.adj = 0)
  legend(min(gene_coords[,1])-55,max(gene_coords[,2]) - 25,legend = names(name2color),pch = 21,pt.lwd = 1,cex=3,pt.bg=as.character(name2color),box.lwd = 0,bg=NA,pt.cex = 4,title = "Cell type",title.adj = 0)
  dev.off()
  
  write.table(ligands, row.names = F, quote = F, col.names = F, file = "results/expressed_ligands.txt")
  write.table(receptors, row.names = F, quote = F, col.names = F, file = "results/expressed_receptors.txt")
  return(list(gene_coords,rec_lig_3,ligands,receptors,wmax,dus))
}

tgconfig::override_params("annotations/lung_params.yaml","metacell")
l = lig_rec_map(mat_id = id,mc_id = id,fig_nm = "results/figure3/interaction_map.png")
gene_coords = l[[1]]
rec_lig_3 = l[[2]]
ligands = l[[3]]
receptors = l[[4]]
wmax = l[[5]]
dus = l[[6]]
save(file = "saved_work/lig_rec_vars.Rda",gene_coords,rec_lig_3,ligands,receptors,wmax,dus)

The ligand-receptor map of lung development pooled across all time points. Genes (ligands and receptors) were projected on a 2D map based on their correlation structure

Immune vs Non-immune interaction map

tgconfig::override_params("annotations/lung_params.yaml","metacell")
mat = scdb_mat(id)
load("saved_work/lig_rec_vars.Rda")
genes = union(ligands,receptors)
m = t(apply(dus[genes,], 1, tapply, as.vector(mat@cell_metadata[colnames(dus), "sorting.scheme"]), sum))
sizes = table(as.vector(mat@cell_metadata[colnames(dus), "sorting.scheme"]))
m = sweep(m,2,as.vector(sizes),"/") * min(sizes)
z = (m[,2] + 10) / (m[,1] + 10)
si_class = ifelse(abs(log2(z)) > 1, ifelse(log2(z) > 0, "immune", "stroma"), "both") # annotate each gene with its identified cell type (immune vs stroma)
si_col = c("gray", "green3", "red2")
png("results/figureS3/si_specificity.png", height=1000, width=1000)
par(mar=c(5,5,5,5))
lim = c(log2(10), max(c(log2(m[,1] + 10), log2(m[,2] + 10))))
plot(log2(m[,1] + 10), log2(m[,2] + 10), pch = 21, cex = 2, bg = si_col[ as.numeric(factor(si_class))], 
    axes = F, xlab = "Non immune", ylab = "Immune", xlim = lim, ylim = lim, cex.lab =2)
abline(coef = c(1,1), lty=2, lwd = 2); abline(coef = c(-1,1), lty=2, lwd=2)
axis(1,cex.axis = 2); axis(2,cex.axis = 2)
dev.off()

## png 
##   2

png("results/figure3/interaction_si.png", height = 1500, width = 1500)
par(mar=c(3,20,3,0),xpd=TRUE)
plot(gene_coords[,1], gene_coords[,2], type = "n", axes = F, xlab = "", ylab = "")
with(rec_lig_3[ rec_lig_3$ligand %in% genes & rec_lig_3$receptor %in% genes & rec_lig_3$interaction == 1,],
    segments(gene_coords[ligand,1], gene_coords[ligand,2], gene_coords[receptor,1], gene_coords[receptor,2], lwd = 1.5,
        col = "gray60"))
points(gene_coords[,1], gene_coords[,2], pch = 21, cex = 4, lwd = 4,
    bg  = ifelse(rownames(gene_coords) %in% ligands, si_col[as.numeric(factor(si_class))], "white"),
        col = ifelse(rownames(gene_coords) %in% ligands, "black", si_col[as.numeric(factor(si_class))]))
legend(min(gene_coords[,1])-85,max(gene_coords[,2]) + 10,legend = c("ligand","receptor"), pch=21,pt.lwd = 1,cex = 3,pt.bg=c("gray","white"),box.lwd = 0,bg=NA,title = "LR type",title.adj = 0)
legend(min(gene_coords[,1])-85,max(gene_coords[,2]) - 50,legend = c("Non-immune","Immune","Non-specific"),pch = 21,pt.lwd = 1,cex=3,pt.bg=c("red2","green3","gray"),box.lwd = 0,bg=NA,title = "Niche specificity",title.adj = 0)
dev.off()

## png 
##   2

Projection of genes activated in the immune (green) and non-immune (red) compartments.
Full and empty circles represent ligands and receptors, respectively. Gray circles represent ligand/receptors non-specific to one compartment.

Differential expression of 604 LR genes between the non-immune (red, x axis) and immune (green, y axis) compartments.
Compartment specificity is determined by two-fold change threshold. LR which are not specific for immune or stromal compartment are marked in gray circles.

Compartment specificity

To identify important cellular communication hubs involved in a large number of interactions between and within compartments, we examined LR expression patterns across different cell types.
In this chunk we calculate the fold change of each gene for each cell type, if it is higher than 2 we include the gene in the interaction map of the cell type.

library(reshape2)
source(src_functions)
load("saved_work/lig_rec_vars.Rda")
genes = union(ligands,receptors)
disp_genes = c("Met","Robo2", "Wnt5a","Sdc2","Tgfb3","Il33","Areg", "Sftpa1","Icam1","Sdc1","Cgn","Cd4","Tnfsf11","Cd247","Il17f","Chad","Il22","Csf2","Il18r1","Il2ra","Il18rap","Il1rl1","Csf2ra", "Csf2rb","Ccl6","Ccl9","Il6","C3ar1","Ccl4","Il6","Il4","Csf1","Hgf","L1cam","Il13","Ccl24","C1qb","Ccl7","Pf4","C3ar1","Sirpa","Ccl3","Csf1r","Tnf","C3ar1","Il6ra","Csf2rb","Spp1")
lin_ord = c("Endothel","Fibro","Matrix","Smooth","Pericytes","Epithel","AT1","AT2","Club","Ciliated","MacI","MacII","MacIII","Mon","Neut", "Baso", "Mast","DC","B", "T", "ILC","NK")
pop_fc = matrix(NA, nrow = length(genes), ncol = length(lin_ord), dimnames = list(genes, lin_ord))
mc = scdb_mc(id)
color_key = unique(mc@color_key[,c("group","color")])
name2color = color_key$color
names(name2color) = color_key$group
color2name = color_key$group
names(color2name) = color_key$color
short_list = rec_lig_3[ rec_lig_3$ligand %in% genes & rec_lig_3$receptor %in% genes & rec_lig_3$interaction == 1,]
for (pop in lin_ord) {
    g1 = intersect(colnames(dus),names(mc@mc[which(as.character(color2name[mc@colors[mc@mc]]) == pop)]))
    g2 = setdiff( intersect(colnames(dus),names(mc@mc)), g1)
    x = rowSums(dus[genes,g2]) / length(g2) * min(length(g1), length(g2))
    y = rowSums(dus[genes,g1]) / length(g1) * min(length(g1), length(g2))
    z = log2((10 + y) / (10 + x))
    spec_genes = names(which(z > 1))
    pop_fc[,pop] = z
    all_int = rec_lig_3[ (rec_lig_3$ligand %in% spec_genes | rec_lig_3$receptor %in% spec_genes) & rec_lig_3$interaction == 1, 1:2]
    all_genes = intersect(rownames(gene_coords), as.vector(as.matrix(all_int)))
    png(paste0("results/figure3/pops/", pop, ".png"), height = 1200, width = 1200)
        plot(gene_coords[,1], gene_coords[,2], type = "n", axes = F, xlab = "", ylab = "")
        points(gene_coords[,1], gene_coords[,2], pch = 20, cex = 2, lwd = 2.5, col = "gray80")
    with(short_list[ short_list$ligand %in% spec_genes | short_list$receptor %in% spec_genes,],
                segments(gene_coords[ligand,1], gene_coords[ligand,2], gene_coords[receptor,1], gene_coords[receptor,2], lwd = 1.5, col = "gray60"))
        points(gene_coords[all_genes,1], gene_coords[all_genes,2], pch = 21, cex = 6, lwd = 6,
                bg  = ifelse(all_genes %in% ligands, ifelse(all_genes %in% spec_genes, name2color[pop], "gray60"), "white"),
                col = ifelse(all_genes %in% ligands, "black", ifelse(all_genes %in% spec_genes, name2color[pop], "gray40")))
    dev.off()

    sub_genes = intersect(spec_genes,disp_genes)
    print(paste(pop,max(z)))
    png(paste0("results/figure3/pops/", pop, "_text.png"), height = 1000, width = 1000)
        plot(gene_coords[,1], gene_coords[,2], type = "n", axes = F, xlab = "", ylab = "")
        with(short_list[ short_list$ligand %in% spec_genes | short_list$receptor %in% spec_genes,],
                segments(gene_coords[ligand,1], gene_coords[ligand,2], gene_coords[receptor,1], gene_coords[receptor,2], lwd = 1.5, col = "gray60"))
    points(gene_coords[all_genes,1], gene_coords[all_genes,2], pch = 21, cex = 6, lwd = 4,
                bg  = ifelse(all_genes %in% ligands, ifelse(all_genes %in% spec_genes, name2color[pop], "gray60"), "white"),
                col = ifelse(all_genes %in% ligands, "gray20", ifelse(all_genes %in% spec_genes, name2color[pop], "gray40")))
        if (length(sub_genes) > 0) {
        text(gene_coords[sub_genes,1], gene_coords[sub_genes,2], sub_genes, cex = 3, col = "black",adj = c(0.5,1.2))
    }
  dev.off()
}

## [1] "Endothel 5.68366839382905"

## [1] "Fibro 4.62778283694834"

## [1] "Matrix 4.621435842016"

## [1] "Smooth 5.30596106147188"

## [1] "Pericytes 3.17098079101486"

## [1] "Epithel 3.01153319478587"

## [1] "AT1 5.11735373586667"

## [1] "AT2 6.45774951565094"

## [1] "Club 6.6287500421409"

## [1] "Ciliated 2.17820180930405"

## [1] "MacI 5.39521617381261"

## [1] "MacII 3.46276249483392"

## [1] "MacIII 4.20292951572517"

## [1] "Mon 3.00068754158198"

## [1] "Neut 6.69675505695389"

## [1] "Baso 5.31178386964661"

## [1] "Mast 4.61201260494894"

## [1] "DC 4.10693058339278"

## [1] "B 4.69537570008458"

## [1] "T 4.19341329149704"

## [1] "ILC 2.63035337822058"

## [1] "NK 6.08251400569485"

# Generating a table of all found interactions between lung cell types
tps = factor(mat@cell_metadata[names(mc@mc),"kinetic.group"], levels = c("E12.5", "E16.5", "E_late", "P_early", "P_mid", "P_d2", "P_d7"))
names(tps) = names(mc@mc)
comb = paste0(color2name[mc@colors[as.integer(mc@mc)]], "@", tps)
names(comb) = names(mc@mc)
m = sc_to_bulk(mat_id = id, comb = comb, choose_genes=F, min_comb = 9)
pops = unlist(lapply(sapply(colnames(m), strsplit, "@"), "[[", 1))
pop_tps = unlist(lapply(sapply(colnames(m), strsplit, "@"), "[[", 2))
short_list = short_list[,1:3]
pop_melt = melt(pop_fc)
pop_melt = pop_melt[ pop_melt$value > 1,]
short_list_2 = merge(short_list, pop_melt, by.x="ligand", by.y = "Var1")
short_list_3 = merge(short_list_2, pop_melt, by.x="receptor", by.y = "Var1")
colnames(short_list_3) = c("receptor", "ligand", "interaction", "transmitter", "transmitter_fc", "responder", "responder_fc")
short_list_3 = short_list_3[,c(2,1,4,6,5,7)]
genes = union(short_list_3$ligand, short_list_3$receptor)
m = m[ genes,]
mm = melt(m)
mm$pop = unlist(lapply(sapply(as.vector(mm$Var2), strsplit, "@"), "[[", 1))
mm$tp = unlist(lapply(sapply(as.vector(mm$Var2), strsplit, "@"), "[[", 2))
mm$id = paste0(mm$Var1, "@", mm$pop)
mm2 = mm[,c(6,5,3)]
m2 = dcast(mm2, id ~ factor(tp, levels = levels(tps)))
m3 = as.data.frame(m2[,-1]); rownames(m3) = m2[,1]

short_list_3$lid = paste0(short_list_3$ligand, "@", short_list_3$transmitter)
short_list_3$rid = paste0(short_list_3$receptor, "@", short_list_3$responder)
short_list_4 = merge(short_list_3, m3, by.x = "lid", by.y = 0)
colnames(short_list_4)[9:15] = paste0("ligand-", colnames(short_list_4)[9:15])
short_list_5 = merge(short_list_4, m3, by.x = "rid", by.y = 0)
colnames(short_list_5)[16:22] = paste0("receptor-", colnames(short_list_5)[16:22])
short_list_5 = short_list_5[,-(1:2)]
write.table(short_list_5, sep = "\t", row.names = F, quote = F, file = "results/figureS3/all_interctions.txt")

ILC interaction map

AT2 interaction map

MAC III interaction map

Basophils interaction map

ligand-receptor pairs

tgconfig::override_params("annotations/lung_params.yaml","metacell")
library(scales)
lig_rec_fn = 'annotations/ligand_receptor_mouse.csv'
proper=function(x) paste0(toupper(substr(x, 1, 1)), tolower(substring(x, 2)))
rec_lig = read.delim(lig_rec_fn, sep = ",", row.names = 1)
ligand_receptor = as.matrix(table(rec_lig$ligand,rec_lig$receptor))
rownames(ligand_receptor) = unlist(lapply(rownames(ligand_receptor),function(x) proper(x)))
colnames(ligand_receptor) = unlist(lapply(colnames(ligand_receptor),function(x) proper(x)))

mc = scdb_mc(id)
mc2d = scdb_mc2d(id)
mat = scdb_mat(id)
color_key = unique(mc@color_key[,c("group","color")])
name2color = color_key$color
names(name2color) = color_key$group
color2name = color_key$group
names(color2name) = color_key$color
lin_ord = c("Endothel","Fibro","Matrix","Smooth","Pericytes","Epithel","AT1","AT2","Club","Ciliated","MacI","MacII","MacIII","Mon","Neut", "Baso", "Mast","DC","B", "T", "ILC","NK")
lfp = log2(mc@mc_fp)
sigs = intersect(intersect(rownames(lfp),rownames(dus)), union(rownames(ligand_receptor), colnames(ligand_receptor))) # list of relevant ligand receptor genes
ligands = intersect(sigs, rownames(ligand_receptor)); receptors = setdiff(sigs, ligands)
rec_lig_2 = melt(ligand_receptor)
colnames(rec_lig_2) = c("ligand", "receptor", "interaction")
rec_lig_2[,1] = as.vector(rec_lig_2[,1]); rec_lig_2[,2] = as.vector(rec_lig_2[,2])
rec_lig_3 = rec_lig_2[ rec_lig_2[,1] %in% ligands & rec_lig_2[,2] %in% receptors & rec_lig_2[,3] == 1,c(1,2)]
rownames(rec_lig_3) = paste0(rec_lig_3$ligand,"_",rec_lig_3$receptor)
genes = union(rec_lig_3$ligand,rec_lig_3$receptor)
mc_mc = mc@mc[intersect(colnames(dus),names(mc@mc))]
m = t(apply(dus[genes,], 1, tapply, color2name[mc@colors[mc_mc]], sum))
sizes = table(color2name[mc@colors[mc_mc]])
m = sweep(m,2,as.vector(sizes),"/") * min(sizes)
q = 9

pairs = c("Il33_Il1rl1","Csf1_Csf1r")
seq10 = seq(0,1,length.out = 10)
gb_mat = outer(seq10, seq10, function(x,y) sqrt(x^2 + y ^2))
rownames(gb_mat) = seq10; colnames(gb_mat) = seq10
gb_mat = melt(gb_mat)
colnames(gb_mat)  = c("green", "blue", "gb_norm")
gb_mat$col = ifelse(gb_mat$gb_norm >0.1, hsv(h = ifelse(gb_mat$gb_norm > 0, 0.25 * gb_mat$green / gb_mat$gb_norm, 0), s = pmax(gb_mat$green,gb_mat$blue), v = 1 - 0.5*gb_mat$green - 0.5*gb_mat$blue), "white")

for(pair in pairs){
  lig = strsplit(pair,split = "_")[[1]][1]
  rec = strsplit(pair,split = "_")[[1]][2]
  a = mat@mat[lig,]
  b = mat@mat[rec,]
  norm_a = rep(0, length(a)); names(norm_a) = names(a)
  norm_a[ a != 0] = as.numeric(cut(a[ a != 0], unique(quantile(a[ a != 0], (0:q)/q)), include.lowest = T)) + 1;
  norm_b = rep(0, length(b)); names(norm_b) = names(b)
  norm_b[ b != 0] = as.numeric(cut(b[ b != 0], unique(quantile(b[ b != 0], (0:q)/q)), include.lowest = T)) + 1;
  green = norm_a / max(norm_a); blue = norm_b / max(norm_b)
  gb_norm = sqrt(green^2 + blue^2)
  cols = ifelse(gb_norm >0.1, hsv(h = ifelse(gb_norm > 0, 0.25 * green / gb_norm, 0), s = pmax(green,blue), v = 1 - 0.5*green - 0.5*blue), "white")
  exp_cells = names(which(a > 0 | b > 0))
  pt_cex = 1 + 0.4 * round((pmax(norm_a, norm_b) - 1) / max(norm_a) * 5)
  png(paste0("results/figure3/pairs/", lig, "-", rec, ".png"), height = 1000, width = 1000)
  layout(matrix(c(1,2,3,4),nrow =2,ncol=2,byrow = TRUE),width = c(200,800,300,700),height=c(700,300))
  top_left_marg=c(5,3,40,1)
  par(mar=top_left_marg)
  image(matrix(1:100, nrow = 10), col = gb_mat$col,axes=F)
  arrows(-0.1,-0.1,-0.1,0.5,xpd=TRUE)
  arrows(-0.1,-0.1,0.5,-0.1,xpd=TRUE)
  mtext("Receptor",side=2,adj=0,line = 1)
  mtext("Ligand",side=1,adj=0,line = 1)
  top_right_marg =c(0,1,1,1)
  par(mar=top_right_marg)
  # plot(mc2d@sc_x, mc2d@sc_y, pch = 8, col = alpha(mc@colors[mc@mc[names(mc2d@sc_x)]],alpha = 0.4), cex=1, axes = F, xlab = "", ylab = "")
  plot(mc2d@sc_x, mc2d@sc_y, pch = 8, col = "gray80", cex=1, axes = F, xlab = "", ylab = "")
  points(mc2d@sc_x[exp_cells], mc2d@sc_y[exp_cells], cex = pt_cex[exp_cells],pch = 21,bg = cols[exp_cells],lwd = 1.5)
  legend_marg = c(1,1,1,1)
  par(mar=legend_marg, font = 2)
  plot(NULL ,xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
  legend("center",legend=names(name2color[lin_ord]),title = "Cell types colors", text.font = 1, pch=19, cex=1.5,col=as.character(name2color[lin_ord]),ncol = 2, bty='n')
  bottom_marg = c(3,2,3,1)
  par(mar=bottom_marg,lwd=3,cex = 1.2)
  ylim = c(floor(-max(m[rec,])), ceiling(max(m[lig,])))
  barplot(m[lig,lin_ord ], col = name2color[lin_ord], names.arg="", ylab = "", ylim = ylim, axes = F)
  barplot(-m[rec,lin_ord],  col = name2color[lin_ord], names.arg="", add = T, axes = F)
  axis(2, at = c(ylim[1], 0, ylim[2]),las =2,font=2)
  mtext(lig,side=3,adj = 0.5,line = 1.2,col = "#408000",cex = 2,font=2)
  mtext(rec,side=1,adj = 0.5,line= 1.2, col = "#800000",cex = 2,font=2)
    dev.off()
}

Dual projection of the ligand Csf1 (green) and its unique receptor Csf1r (red) on the single cell map from Figure 1.
Colors indicate expression quantiles. Bar plots indicate ligand and receptor normalized expression per 1,000 UMI across cell types.